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Numerical Solution of the Navier-Stokes Equations 
for Blunt Nosed Bodies in Supersonic Flows 

By 

Z, U. A. Warsi, Krishna Devarayalu and J. F. Thompson 
Department of Aerophysics and Aerospace En g in eering 
Mississippi State University 
Mississippi State, MS 39762 

Abstract 

A time dependent, two-dimensional Navier-Stokes code employing the 
method of body-fitted coordinate technique has been developed for super- 
sonic flows past blunt bodies of arbitrary shapes. The bow shock ahead 
of the body is obtained as part of the solution, viz., by "shock capturing." 
A first attempt at mesh refinement in the shock region has been made by 
using the forcing function in the coordinate generating equations as a linear 
function of the density gradients. This technique displaces a few lines 
from the neighboring region into the shock region. Numerical calculations 
for Mach numbers -2 and 4.6 and Reynolds numbers from 320 to 10^ have been 
performed for a circular cylinder with and without a fairing. In this 
report results of Mach number 4.6 and Reynolds number 10*^ for an isothermal 
wall temperature of 556“K have been presented in detail. 


iil 





List of Symbols 


A 



E. 




f 


F 

s 

G 

G. . 

ij 

H 


= CTu, Eq. (3.16) 

= Amplitude Factor, Eq. (2.12) 

= Amplitude Factor, Eq. (2.11) 

= Vector, Eq. (3.4) 

= ov, Eq. (3.16) 

= Amplitude Factor, Eq. (2.12) 

= Amplitude Factor, Eq. (2.11) 

= Specific Heat at Constant Pressure (Dimensional) 

= Pressure Coefficient, Eq. (3.29) 

= Deformation Tensor, Eq. (3.4) 

= Decay Factor, Eq. (2.12) 

= Decay Factor, Eq. (2.11) 

/ 

= ^ i = 1,2 , j = 1,2 , Eq. (3.10) 

= Non-dimensional Specific Internal Energy 
= , Eq. ,(3.10) 

= Decay Factor, Eq. (2.12) 

= pecay Factor, Eq. (2.11) 

= Scalar Function (Surrogate) 

= Numerical Vector, Eq. (4.11) 

= Det (g..)» Eq. (A-4) 

= Covariant Components of the Metric Tensor, Eq. (A-1) 

= Contravariant Components of the Metric Tensor, Eq. (A-2) 
= Numerical Vector, Eq. (4.11) 

- 8y//i 

= a(u)^, Eq. (3,16) 


iv 



j J • 

i 


I 

n-IAX, JIIAX 

I T 
Bl* B2 

^UC’ ^LG 
J 

k*- 

•K 


M 


M 


00 


n 

n,n’ 


n' 

P 

P 

P 

P 

r 

Q 

Q 

R 


“ Used in Index Notation 

= Idem Tensor; Equal to when Referred to Cartesian Coordinates 
and or gij when Referred to . Curvilinear Coordinates, 

Eq. (3.4) 

=* Integers 

= Maximum Values' of I and J Used 

= Values of I where the Cut Meets the Body, Fig. 5 

= Integers Representing I Values on the Upper and Lower 
Parts of the Cut Respectively, Fig. 5 

= Integers; also Jacobian where Stated 

= Conductivity (Dimensional) 

= Xdivv, Eq. (3.4) 

= Number of Points Involved in Attraction Along the ^-Direction, 
Eqs. (2.11), (2.12) 

= cr(v)2, Eq. (3.16) 

= Local Mach Number, Eq. (3.23) 

= Freestream Mach Number 

= Unit Normal Vector, Eq. (B-21) 

= Number of q-Lines Involved"in n-Line Attraction, Eqs. (2.11), 

( 2 . 12 ) 

= auv, Eq. (3.16) 

= Nondimensional Pressure 
= /g p, Eq. (3.10) 

= Function Appearing in Eq. (2.9) 

)i*C* 

= ^ , Prandtl Number, Eq. (3,4) 

k* 

= (E+P)u, Eq. (3.16) 

= Function Appearing in'Eq. (2.9) 

= (E+P)v, Eq. (3.16) 


V 



R, 




R 

e 

R* 

n. 

sgn 

t 

T 

u 

V 

U,V 

i 


V. 

X 


i 

V . 

>3 

w 

(x,y) 

a 

o 

Jk 

Y 

C,n 


^k*\ 

0,^1, i|; 

'V 

P 


= Numerical Vector, Eq. (4.16) 


2P*V* R* 
“ oo n 


, freestream Reynolds Number 


= Nose Radius (Dimensional) 

= Sutherland’s Constant, Eq. (3.7) 
= Sign Function, Eq. (2.11) 

= Unit Tangent Vector, Eq. (B-19) 


_ ^ Nondimens ional Temperature 

= v^, Contravariant Velocity Component 
=. v^, Contravariant Velocity Component 
= Cartesian Velocity Components 
= Contravariant Components 
= Covariant Components 
= Y /V^, Nondimensional Velocity Vector 


= Covariant Derivative 
= Numerical Vector, Eq. (4. 11 ) 

= Nondimensional Cartesian Coordinates 
= Cartesian Coordinates (x^^ = x, = y) 
= Nondimensional Constant, Eq. (3.4) 

= Christoff el Symbols, Eq. (A-7) 

= c*/c^. Adiabatic Exponent 
= Transformed Coordinate System 
= Any Curvilinear Coordinate System, ('|^ 


= n) 


= Values of 5 and n for Varying k, Eq. (2.11) 
= Functions, Eqs. (A-16) - (A-18) 

= Total Energy, Eq. (3.4) 

j ^ 

= p /p", Nondimensional Density 


Vi 



a 


o 

y 

X 


vl 

T 

e 

u 




6 .. 

13 


= V^p, Eq. (3.10) 

= Stress Tensor, Eq. (3.4) 

- y /\x^y Nondimensional Viscosity 

= X*/ii*, Nondimensional Bulk Viscosity; Also a Surrogate 
Variable where Stated 

= Functions, Eq. (3.14), (A-24), (A-25) 

= Tensor, Eq. (3.4) 

= 1/R 

e 

= Vorticity; Also as Acceleration Parameter where Stated 
= Kronecker-6 


Operators ; 

div = Divergence 

grad = Gradient 

( )"^ = Transpose 


Subscripts ! 

= Denotes Vector 

5,Ti = Partial Differentiation 

■» = Frees tream Value 

w = Wall Value 

ij - .Covariant Tensor 

= Covariant Derivative 

i = Covariant Vector; when used V7ith x or then no Tensorial 

Significance 

s = Stagnation Value 

Superscripts : 

"Vj = Denotes a Tensor 

* = Denotes a Dimensional Quantity 

ij = Contravariant Tensor 

i = Contravariant Vector 


vii 



!• Introduction 


Niamerical solution of the full compressible Navier-Stokes equations for 
a blunt body placed in a uniform supersonic stream is of much theoretical 
and practical interest. The problems of viscous compressible flow have been 
tractable only in the last decade because of the continuing advances in 
the computer facilities as well as in the numerical methods. A detailed 
summary of various numerical methods which have been used in the calculation 
of viscous compressible flow problems has recently been given by Peyret 
and Viviand [1]. A recent review on blunt body problems by Rusanov [2] 
is also of interest. 

It has been established by experiments that in front of a blunt body 
placed in a supersonic flow, there appears a detached bow shock wave which 
separates the flow disturbed by the body from the undisturbed flow. Behind 
the shock wave there is a subsonic region bounded by it, the body, and the 
sonic lines emanating near the body. The region behind the sonic lines is 
again supersonic with the exception of interior shocks and near-body wake. 
The existence of the subsonic region, with the no-slip condition to be 
satisfied at the body surface, requires the consideration of full viscous 
equations. The subsonic region also establishes a direct relationship 
between the shapes of the nose part of the body and that of shock wave 
ahead of it. 

The present investigation is concerned with the numerical solutions 
of the complete Havier-Stokes equations for two-dimensional blunt arbitrary 
shaped bodies placed in moderately high freestream Mach number flows. One 
of the novel features of the present research is the use of body-fitted 
numerically generated coordinate systems as developed by Thompson, et. al. 
[3] and Thames [4]. The coordinate system used in this paper (Ref. to 
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Figs. 1, 2, and 3) establishes beyond doubt the versatility of the method 
and the capability of adjustments of the coordinates in regions of interest, 
such as in the vicinity of a shock or a body surface. Further, the overall 
problem of analytical development has been rendered more general by using 
the method of tensor analysis. Thus the extension of the method to three 
dimensions is now purely formal. 

The body-fitted coordinate generation method allows the coordinate 
lines to be coincident with all the boundaries of a general multiply- 
connected region including the boundaries formed by solid walls and the 
external boundaries. Thus with this procedure the numerical solution of 
the Navier— Stokes equations may be obtained on a fixed rectangular field 
in the transformed plane without any specification of the mesh sizes, 
cf. [5]. Further, no interpolation of the flow variables is required 
regardless of the shape of the physical boundaries or the spacing of the 
curvilinear coordinate lines in the physical plane. 

As noted earlier, the flow in the subsonic region is governed by the 
shapes of the nose part of the body and that of the shock wave ahead of it. 

By using the body-fitted coordinate system, the problem of exact specification 
of the boundary conditions on the nose part of the body is completely 
resolved. In Sections 2 and 4, both the theory and numerical generation 
of coordinates are discussed. 

The computational domain is now limited upstream by a boundary located 
at a certain distance ahead of the expected appearance of the bow shock in 
the steady state. The outer boundary is taken to be a hyperbola placed 
at a distance about 5.0 of the expected stand-off distance calculated by 
Van Dyke [6] and. Billig [7]. The downstream boundary is a circular arc 
of radius equal to 5 times the nose radius of the body. An implicit finite- 
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difference approximation for the non-steady Navier-Stokes equations written 
in curvilinear coordinates is used to advance the solution from an arbitrary 
specified solution at time t = 0, Since the position or shape of the shock 
has not been prescribed priori, the shock transition region is obtained 
as part of the solution, i.e,, by "shock capturing," as the flow develops. 
Section 3 describes all the pertinent equations used In the flow calculations 
while Section 4.2 describes the solution algorithm. 

As is usual with any finite-difference solution of the compressible 
ITavier- Stokes equations, the problem of nonlinear oscillations, particularly 
when a shock exists, is a dominant one. These oscillations or wiggles must 
be damped by using dissipative finite-difference schemes. McCormack [8], 
Boris and Book [9], Vliegenthart [10] and others have developed various 
techniques to damp out the nonlinear oscillations. In the present research 
we have used both the FCT routine of Ref. [9] and the Shuman filtering of 
Ref. [10]. It has been found that Shuman filtering works quite well for 
all the cases considered in this paper, though it has the tendency of 
lessening the high gradients in the shock region. - This aspect of the problem 
is discussed in Section 4.2.3. 

A technique to concentrate the coordinate lines in the vicinity of 
the shock has been developed. When a quasi steady-state solution has been 
obtained, the density gradients already known across the shock are then used 
to re- generate the coordinates and the Navier— Stokes solution is performed 
on the new coordinates. Though this scheme works quite well, however, 
because of the limited storage capacity, not many lines could be displaced 
in the transition region. This aspect is discussed in Section 4.2.5. of this 
paper. 
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A series of results on a circular cylinder and a circular cylinder 
with aft-body have been obtained for several Reynolds and Mach numbers. 

The wall pressure ratio on a circular cylinder as reported by Tannehlll and 
Holst [11] has been compared with the present solution and the two match 
fairly well. All numerical results are discussed in Section 5 of 
this paper. In all cases considered the Knudsen number, Air/2 M /R , is 
sufficiently small to satisfy the continuum hypothesis. 
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2, Numerical Generation of Curvilinear Coordinates 


The accuracy of numerical solutions of the partial differential equations 
of mathematical physics depends very strongly on our ability to impose the 
requisite boundary conditions as accurately as possible* For problems of 
fluid flows past finite bodies > since the boundary values are prescribed 
on a closed curve forming the body contour (in three dimensions on a closed 
surface) , and also at another closed curve enclosing the given body fomning 
the outer or external boundary, it is natural to envisage a coordinate 
system in which one coordinate passes exactly through the body contour, 
and the other passes through the external boundary contour. The idea of 
numerical coordinate generation is to- fill the region with intersecting 
coordinate lines enclosed by the body and the external boundary contours in 
the physical (x,y) or (x,y,z) space. 

The preceding ideas in one form or another have been used by Winslow 
[12], Barfield [13], Chu [14], Amsden and Hirt [15], and Godunov and 
Prokopov [16]. However, the whole concept has been used in a much organized 
manner to provide a number of solutions in fluid mechanics by Thompson, et. 
al*, [3], [4], [5], [17I, [18], [19], [20]. 

Let 5 ” €(x,y) and n - n(x,y) be two continuously differentiable 

functions of the Cartesian coordinates (x,y). Further, let n = n = constant 

0 

be the body contour, while n == = constant be the external boundary con- 

tour, The region < n < n^must now be filled by intersecting, coordinate 
curves ^ = constant and n = constant. Becuase of the closed region under 
consideration it is natural to specify the determining differential equations 
for 5 and n es elliptic equations to be solved under the proper boundary 
conditions for g and n at the body and at the external boundary. Since the 
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simplest elliptic equation is the 'Laplace equation, we then pose the problem 
of solving the Laplace equations for ^ and n with x and y as independent 


variables under the Dirichlet boundary conditions. Let be the curve 
defining the body contour ^ ~ a-nd T2 he the curve defining the outer 
boundary ri = in the xy-plane as shown in Fig* 1* The elliptic boundary 
value problem is then 



< 

ro 

11 

0 

( 2 . 1 ) 


<j 

10 

11 

0 

( 2 . 2 ) 

on T^: C 

=,f^(x,y), n = 

(2.3) 

on r^: S 

= f„(x,y), n = 

(2.4) 


The solutions of Eqs. (2*1) and (2*2) under the boundary conditions 
(2.3) and (2*4) can conveniently be obtained in those cases when q and n 

O' “ 

can be specif ied' by simple analytic methods (such as a circle, ellipse, 
etc*). To obtain coordinates for arbitrary shaped bodies, it is convenient 
to transform the Eqs. (2.1) and (2.2) such that x and y are the dependent 
while 5 and n the independent variables. This transformation can easily 
be performed for either two or three-dimensional coordinates by the method 
of tensor analysis and is detailed in Appendix B* Referring to Eqs* (B-13) 
and (B-14) , we find that Eqs. (2.1) and (2.2) are equivalent to 


'22 " 

2gj^2 

*5n + ®11 - 0 

(2.5) 

22 “ 


=^5n Su/nn ° “ 

(2.6) 


where the variable subscripts denote partial differentiations and the g 

Ij 

is the metric tensor defined in Appendices A and B. The boundary conditions 
are now 


On r*; 


y - 

(2. 

.7) 

On r*: 


y.- Gj - 

(2. 

.8) 
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where as shown in Fig. 2, F*- and F* are the images of the body and the external 
boundary contours in the ^Ti-plane, 

The geometrical meaning of the transformed equations (2.5) and (2.6) is 
that the body and the external body contours in the xy-plane have been mapped 
on to the ^n~plane which is rectangular. In other words, we can say, 
that the contours in the xy-plane have been opened up to form the straight 
lines n = = constant, and H = constant in the ^n-plane-. This can 

be achieved by imagining a cut connecting the body and the external boundary 
in* the xy-plane as shoxra in Fig. 1. , such that all functions and their 
derivatives are continuous in crossing the cut. Since a- cut line is a part 
of the field, no boundary conditions can be imposed on F* and F* of Fig. 2. 

the appearance of n and t| in Eqs. (2.7) and (2.8) is now purely 
symbolic, denoting the names of the body and of the external boundary 
respectively. Given the body and the external boundary contours, we can 
always establish the values of x and y either graphically or analytically 
for any desired distribution of ^"Values. The q-values can be chosen 
arbitrarily to form rectangular meshes in_the Sq-plane. 

Equations (2,5) and (2.6) are the basic equations for the generation 
of coordinates. To have a control over the spacing of the C aud n lines, 
we envisage another general transformation, say from to C^'and U'*. 

Retaining the 5,ri notation, the equations take the form 

*22 ■ ^*12 + *11 *.in % 

*22 ^55 - 2*12 *11 ^nn ’ % * 

For details on the above derivations refer to [ 51i, The same form of 

equations are obtained if one starts considering the Poisson equations in 

place of Eqs. (2.1) and (2.2) and inverts the transformation from x,y to 

5 , n as independent variables . 


(2.9) ; 

( 2 . 10 ) 
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The function P and Q are to some extent arbitrary and can be chosen 
in various ways to have a desired distribution of coordinates in a given 
region. In the present research we have made p and Q to depend on the 
density gradients to contract the coordinates in the region of the shock. 


The chosen forms of P,Q are [17] 


n^ 


= 8 

E 

k=l 

K 

sgn 

(?- 



m^ 





+ g 

E 

1=1 


sgn 

(?- 



n 





= 8 

E 

k=l 


sgn 

(q- 



exp 

exp (-E^R^) (2.11) 

exp (-Dj^In-n^|) 


where 


• ™ 

+ S sgn (n-ri^) exp (“Ej^R^) 

X — i 

® " ®11 ®22 ~ ^® 12 ^^ ^ 


(2,12) 


(2.13) 


\ (n-Hj^)2]^^^ (2.14) 

The first terms on the right hand sides of both (2.11) and (2.12) are used 
in the line attraction, while the second terms in both equations are used 
for the point -at traction. Various terms which appear in these equations 
have been defined in the "List of Sjnnbols". 

In the present research, we have used only the .point attraction term 
of Eq. (2.12) to concentrate the coordinate lines near the shock. The 
computer program calculates the density difference for all q-lines in the 
region of the shock. The amplitude factor B. thus changes according to the 

Jv 

position (5,n) and is defined as 

B^ = (constant) (p 2 ~Pj^)/p^ (2.15) 
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where the subscripts 1 and 2 denote the respective values in the front and 
behind of that shock which has been computed without coordinate contraction. 

The constant appearing in (2.15) is selected only once by trial and error 
and retains’ the same value for all g and n positions. 

The method of numerical coordinate generation offers much freedom in 
the orientation of both the 5 and r\ coordinates in the physical xy-plane. 

For example, the r\ = const, lines can be chosen to go round the body ,as 
shown in Fig. 1, or they may not be chosen to form a complete circuit as 
shown in Fig. 3. However, a choice has to be made in advance of computing 
the coordinates, because the resulting orientations of the body segment, 
the cut lines, or the re-entrant segments, and the outer boundary segments 
in the 5TV“plane depend on this choice. In the present research we have 
chosen the coordinate orientation as shown in Fig. 3, in which the front 
outer boundary is a hyperbolic arc and the rear outer boundary is a circular 
arc. Figure 4 shows the corresponding segments orientation in the ?n-plane. 
This type of segments orientation requires much care in the computer programing 
for both the coordinate generation and for the numerical solution of the 
Navier— Stokes equations. In Section 4 we have discussed the finite-difference 
approximations of all the equations. 
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3. Formulation of the Problem 

For solving the blunt body problem in a general curvilinear coordinate 
system, we consider the nondimensionalized Navier^ Stokes system of equations 
in the invariant tensor form. The conservation equations are 
Mass Conservation: 


^+div(py) = 0 (3.1) 


Momentum Conservation: 


(py) ’+ div T = 0 

Energy Conservation: 


(3.2) 


where 


t 


— + div b = 0 (3.3) 

T = pyy + pi “ ea 
b = (’f+p)v - e5*v - a^y grad T 
'? = pe + p |yp 
5 = KI + d 

K = X div V 


t vectors and tensors are denoted by using - under and above a symbol 
respectively. 
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d = ydef y = y[grad y + (grad y)^] 
e = 1/R 

e 

“O " 

(3.4) 

The nondimensionalized equations (3.1) — (3.4) have been obtained 

by referring all lengths to the diameter 2R*; velocity vector, density, 

viscosity, temperature, enthalpy and pressure to the freestream values 
* * * * * * * 2 

» ’ Poo ^00 ’^fispectively. Assuming the gas to 

be caloricaily and thermally perfect, the equation of state in the non- 
dimensional variables is 



(3.5) 


where y is the adiabatic exponent. Similarly, the nondimens ional tempera- 
ture is given by 


T = y(y-1)m 2 (^ - I |v|2) 

The relation between temperature and viscosity is provided by the Sutherland 
formula, which in non-dimensional form is 


where 


* 



* 

T 


00 


U 


(l+S^)T^^^ 


* 

= llO^K . 


(3.7) 
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Since 4' is a function of T, hence the system of equations (3.1) - (3.7) 
form a closed system of equations provided X = A(u) is also prescribed. 
In the present formulation we have used the Stokes* condition 


3A + 2p = 0 


(3.8) 


as the required relation. 

The boundary conditions for the system of equations (3.1) - (3.3) 

are 

t I 3T 

at the body surface: |y| = 0, T = T^, or specified 


't'.. = 


T p 
w w 




Py.t 2 


w 


yK 


at freestream infinity: ly| =1, p = 1, T = 1 


(3.9) 


P = 




= 1 + 


2 y(Y“1)M^ 


In (3.9) the subscript w denotes the wall condition. The density is 
not known in advance but must be obtained by the equations themselves. 

Since the governing equations are of parabolic-elliptic type, we 
therefore need to specify, the outflow boundary conditions. In place of 
specifying the derivative conditions, we have used the complete solution 
of the Euler equations to specify the downstream boundary condition with 
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the stipulation that at the outer downstream boundary the effect of viscosity 
is negligible. 

3 . 1 Transformation of the E quations; 

The governing equations (3.1) - (3.3) are now transformed to a general 
coordinate system 5^(i = 1,2) by the usual method of tensor analysis.^ Using 
the summation convention on repeated indices, and introducing new dependent 
variables 


o = p, P = p, d^^*, E = (3.10) 

the equations are 


1^ + (av^) = 0 (3,. 11) 

~ (ov^) + (ov^v^ ) + r.^ av^ 

9t ggj jk 

- p 4gj^ + a/T 

+ + 4 + D^^)} (3.12) 

3g3 JK 

-|| + {(E + P)v^} = e ~~ (/g~ v^) + a- (pv^ -^) (3.13) 

° 3?^ 3?^ 

where is the transformation metric tensor , and 

+ p(g^^g^^ (3.14) 


^ (P/^) 
3r 


iRef er to Appendix ‘ A. 
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In Eqs. (3.10) - (3.14), the superscript indices denote the contravariant 


components, while a comma denotes the convariant differentiation, i. 


i 9V-" , „i k 

V . = r + r ,, V 

,3 3^3 3k 

where is the Christoff el symbol. 

Equations (3.11) - (3.13) are applicable both to two and three 
dimensions. In two dimensions, writing for brevity 

5^ = 5, = n, = u, = V 

A = ou, B = ov, H = o(u)^, M = o(v)^, N = auv, 

Q = (E+P)u, R = (E+P)v 
the equations become of the form 

+ M = 0 




aa 



at 


aH 

+ 3N 

at ^ 

9? 

^ 9n 

M + 

9N 

aM 

•at ^ 

9? 




9E 



at 


.1 f 
11 


I 

11 


'1 N + t 
12^ 22 


■'2 jT + p2 j, 
12“ 22 


ww , vxv , 

The expanded form of various terms appearing in Equations (3.17) 
are given in Appendix A. 

The boundary conditions for Equations (3.17) - (3.20) are 

3T 

at the body surface: u = v= 0, T = T,or (— ) specified 

w on w 


E = 


T o 
w w 


^ Y(Y“l)M:f; 


^w == 


e* 9 

(3.15) 
space 

(3.16) 

(3.17) 

(3.18) 

(3.19) 

X3.20) 

(3.20) 
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at freestream infinity: 


a -> /g 


I* »^/yM^ 

E ^ i^d + ^ 


-) , 


^ Y (y-DM^ 

where a 5 or subscript implies partial differentiation. 


(3 


The local Mach number M is given by 

Xr 




where 


and 


lv|2 = g..(u)2 + + §22^^^^ 


^11’ 


(3 


(3 


T = y(y-1)m 2 (■“ - y 1 y|2) (3 

The relations between the local Cartesian and the local contravariant 
components of the velocity vector v are 

TJ = ux^ + vx 

K n 

V = uy^ + vy^ . (3 


The vorticity m is given by the formula 


_ 1 /’'2 *’^ 1 , 

“ T~) 

3n " 


(3 


where v^ are the covariant components of v, which are related with the 


contravariant components as 


, v^ = g,, u + g, 


’^i ■ ®ij " " Hi “ ‘ H2 


The pressure coefficient is defined as 


C = 
P 


* 

P - 


~ P* V*2 

z 00 CO 


= 2(p - 


(3 


(3 


.22) 

.23) 

.24) 

.25) 

.26) 

.27) 

.28) 

.29) 
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4. Finite Difference Approximation and the Solution Algorithm 

In this section we present the finite difference approximations and the 
numerical methods used both for the solution of the elliptic system, Eqs. 
(2.9) - (2.10), and for the complete Navier-Stokes system, Eqs. (3.17) - 
(3.20). 

Before we proceed on the pertinent method of solution, it is Important 
to mention that in the method of body-fitted coordinates it is superfluous 

r * 

to specify the step sizes A? and An. tf IMAX and JIIAX represent Integers 
for the maximum numbers of C and n points respectively in a field, then 
this input and the desired contraction controlled by the amplitude and decay 
factors of Eqs. (2.11) - (2.12) decide the variable mesh sizes to be obtained 
by solving the generating system (2.9) - (2.10). This aspect has been 
thoroughly discussed in [ 5 ] . Thus the main utility of numerically generated 
body-fitted coordinates actually lies in the availability of meshes or nets 
in the ^n-plane on which the Navier-Stokes equations are to be solved 
without specifying the step sizes. Further, the variations both along the 
?-and n-coordinates , are labeled by the -consecutive integers in the range 
1 < I < IMAX and 1 < J < JMAX. 

4.1 Numerical Generation of Coordinates ; 

The solutions of Eqs. (2.9) - (2.10) have been obtained- by the Gauss- 
Siedel method with -successive over relaxation (SOR) under the prescribed 
boundary values for x and y on the body and the .external boundary contours, 
along with the prescribed values of IMAX and Jt-IAX. The spatial derivatives 
are approximated by the central differences 



(4.1.) 


(4.2) 
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where X is a surrogate variable and I, J denote positions along ? and n 
coordinates respectively. Similar expressions are obtained for X^, X^^ and 



The solutions of Eqs. (2.9) - (2.10) yield x and y for the whole flow 

field as functions of g and q. This data is then used to generate the 

derivatives x_, x , y„, y , the metrical coefficients g.., the determinant 

C’ n S n 

i 

g, and the Christoffel symbols 

It was mentioned in Section 2 that the orientation of the coordinates 


of the form shown in Fig. 3 requires due care in obtaining derivatives on 
the cut line. Referring to the schematic shown in Fig. 5 with the x-axis 
oriented along the cut line, let I and I _ denote the integral values of 

ULi LiO 

I on the upper and lower parts of the cut respectively- Thus 


’^LC + ^OC = ^ • 


(4.3) 


Equation (4.3) establishes the following correspondence between and I. 


LC 


1 

2 


Corresponds to 


UG 


^UC 

IMAX 

IMAX-i 


IMAX-I+1 


B1 B2 (4.4) 

where and represent the same point of the body reached by the lower 
and upper parts of the cut respectively. Obviously 


Ib2 = (IMAX + 1) - 13^ . 
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From (4*4), we conclude that 




(4.5) 


and 


x(I^^+l,J) = x(I^^-l,J) (4.6) 

The first derivatives on. the lower and upper parts of the cut are 

^\hc = 2 

VlC = 2 

where 2 < - 1 . 

(x^)^^ = ~ [x(I^p+l,l) - x(Iy^-l,l)] 

” 2 ~ (4.8) 


where 


"b 2 1 ^ "uc ^ • 


Using (4.6) in (4.8), we find that 




= 0 






(4.9) 
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1 

From the definition of F 

Jh 


given in Eq. 


(B-6) we find that on the cut line 




(4.10) 


4.2 Navier-Stokes Algorithm ; 

The Navier-Stokes equations (3.17) - (3.20) can be put in the numerical 
vector form as 


where 


8w 3F 


3G 


3t 3? 3n 


w 



o 

A 

B 

E 


(4.11) 


(4.12) 



(4.13) 



(4.14) 


H 



H + 2V^2, N + 1^2 M - 6 
r^i H -H 2Vl^ N + r|2 M - 
- 'I' 


(4.15) 
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We now discretize Eq. (4.11) by a fully implicit difference approximation. 
The time derivative is approximated by a first-order backward difference 
at n+1, where n is the time step of size At, while the spatial derivatives 
are approximated by central differences. The system of equations are solved 
by point-SOR, which are 


n+1 _ n+l(p) , 

~I,J ~I,J “~I,J 


(4.16) 


^ the relaxation parameter and the superscript (p) denotes values 
at the previous iteration. The function R is 




•I,J =I,J- 2'I,J 


2 ~I+1,J. ~I-1,J 


* Si,J+l - - <“> Si,j 


where the values of w in G and H are those which are the most recent 
values available from the previous iteration. Fully expanded forms of 
(4.17) are given in [21]. 


4.2.1 n-Derivatives on the Cut ; 

To find the n-derivatives on the cut, we refer to Fig. 6. The point 1 
of the physical plane goes on to the point x of the transformed plane, while 
point 2 remains at its position. Thus in principle, the function value at 
the fictitious point shown under the lower cut must be replaced by the 
function value at the point 2 on the upper cut. Now two cases arise 
depending on whether the function is a scalar (like pressure, temperature, 
etc.,) or a vector. 
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For a scalar f, the first n-derivatives are 






Thus 


(M) = _ (Jil) 

^LC.l ^UC,1 


(4.18) 


Similarly, the second ri-derivatives are 


( 3 ! 1 ) 

3n2 ^nc>l - 

But f is a scalar, so that 


f(lLc.l) = f(I„c,l) 

hence both (4.19) and (4.20) represent the same value. 

To find the n-derivatives of a directional quantity u on the cut, we 
need the value of u at the fictitious point. Since in the physical plane u 
in the lower part of the cut is directed opposite to that on the upper part 
of the cut, hence 




= (— ) 

luc.i (^•21)- 

The same holds for v, x^, x , y_ and y . 

n • ? n 

Based on the preceding analysis it is easy to show that either for a 
scalar function f or a vector function v the derivatives across the cut are 
continuous . 
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4.2*2 Calculation of Wall Density ; 


The wall density is calculated by evaluating each term of the continuity 
equation (3.17) at the body surface. Denoting by J=1 the body surface, we 
have 


"=1,1 


Using a three-point forward difference approximation for the right hand 
side, we obtain 


n+1 n At 

^ 1,1 “ *^ 1,1 “ 2 


(<2 “ ® t 3 > 


(4.22) 


where = 0 . 

Though Eq. (4.22) is fully implicit, nevertheless its use at the 
trailing edge point always produces unrealistic density values. To 
circumvent this difficulty, we used Eq. (4.22) at ail points of the body 
except at the trailing edge point where we have used an explicit scheme 
based on the leap-frog method 


Jia " . _ 


n+1 

iq 


n-1 


2At 




) 


where ‘0’ is a fictitious point. Using a' second order extrapolation in 
space and time, we get 


®I,0 “ “ h,2 


Thus the wall density is obtained by the expression 


A,u /-T.^ I r.n-2,. 

■"l.! “l,! ■ '®I,2 + ®I,2> 


( 4 . 23 ) 
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4.2.3 Nonlinear Instability 

In the calculation of compressible flows, several types of nonlinear 
instabilities are encountered. Among these, the most dominant is due to 
the difference approximation to' the convective derivative. It is a matter 
pf experience that the convective instability can be avoided by introducing 
some dissipation in the difference approximation of the differential equation 
being solved. In this context, McCormack [8] has used a fourth-order 
damping term for his explicit schemes. Boris and Brook [9] have developed, 
a flux-corrected transport (FCT) technique which is quite- efficient for use 
in the continuity equation. Vliegenthart [10] and Harten and Zwas [22] 
have used the "Shuman filtering" to supress the convective instabilities. 

In the present research we have tried both the FCT of Ref. [9] and 
the Shuman filtering of Ref. [10]. Though the Shuman filter adds more 
dissipation than desired, particularly near the shock, it always produced 
wiggles-free solutions for all regions of the flow field. The application 

I 

of Shuman filter amounts to replacing the vector defined by 

j = w“ J J, in (4.17) by. 


“ [ w^ I" w^ *1“ w^ w^ i 4 1 ( 4 . 24 V 

8 I+1,J - I-1,J ~ I,J+1 ~ I.J-1 ^ ^ ~ I,jJ 


4.2.4 Downstream Boundary Conditions ; 

For the solution of the parabolic-elliptic system of equations (3.17) - 
(3.20), beside the boundary conditions at the body surface and -at the upstream 
free boundary (Eqs. (3.22)), it is also necessary to specify a proper set of 

conditions at the downstream boundary. To obtain these conditions, we have 

' \ 

placed the downstream boundary at a sufficient distance where the viscous 
dissipation is negligibly small and have solved the conq)lete Euler ^s equations 
for the whole boundary. The computer program has been structured in such 
a way that is solves both the Navier-Stokes and Euler’s equations simultaneously. 
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4.2.5 Coordinate Contraction Near a Shock: 


The capability of attracting the coordinate lines to other pre- 
designated coordinate lines or grid points exists at present and through 
the application of this technique to blunt body flows with strong shocks 
an effort was made to concentrate in, and define, the region of the shock. 
The magnitude of concentration is controlled by the factors and 
defined in Eq. (2,12). These coordinate control factors were expressed 
as functions of the local density gradients across the shock as shown in 
Eq. (2.15). The equations for the generation of the coordinate system 
(Eqs. (2.9), (2.10)), using coordinate control, were solved, as well as 
the Navier-Stokes equations, when a quasi steady-state has been reached. 

t 

This process was repeated after a pre-assigned number of time steps, 
usually 40. Thus the coordinate system was refined in the region of the 
shock and moved with it. This measure reduced the. need of having a very 
refined mesh in the entire computational domain by providing refinement 
only in the region through which the shock happened to be passing at any 
given time. However, it must be noted that this refinement near the shock 
was achieved at the expense of the accuracy near the wall where coordinate 
contraction is always needed to resolve the boundary layer. 


24 




5. Technique of Numerical Solution 

The first step in the numerical solution of the transformed N-S equations 
is the determination of the computational dor/iain so that the appropriate 
boundary conditions may be prescribed around it* In the present case of 
supersonic flow, the flow field remains unperturbed upstream of the bow shock 
wave, and the computational domain is limited upstreanl by a boundary located 
at a short distance ahead of the bow shock* The "stand-off^* distance of this 
detached shock, for a given freestream Mach number, is estimated using 
empirical equations [ 7 ] . On the upstream boundary the uniform flow con'- 
ditions were used as boundary conditions* Particular care had to be taken 
to ensure that the bow shock did not cut across any segment of this upstream 
boundary. On the downstream boundary the boundary conditions varied with 
time and were determined by solving the Euler equations (cf. Section 4)., 

The computational domain, and the profile of the body in it having been 
determined, the next step was the numerical generation of the coordinate 
system which has already been described. The cartesian coordinates of each 
of the mesh points in the entire computational domain having been determined 
and stored, the coefficients that occur in the Navier- Stokes equations 
due to transforming them into general curvilinear coordinates could now be 
calculated and stored in a file. 

The actual solution of the N-S equations now starts with an assumed 
initial guess of the solution for the entire computational domain. These 
initial conditions need not necessarily be physically realistic and when 
they are not, the transient solution has no physical meaning* In the 
present case the initial conditions chosen for the whole flow field were 
the uniform flow conditions that were prescribed on the upstream boundary* 

It was however found that this could not be done if the freestream Mach 



number was very high, or if the isothermal temperature condition prescribed 
on the body was far different from the freestream value. The finite- 
difference scheme chosen was the S.O.R. which is an implicit scheme. The 
value of the optimum acceleration parameter for all the equations, i.e., 
continuity, momentum and energy, was determined, by trial, to be 0.9. In 
approximating the convective derivatives of these equations, the average of 
the product (A.O.P.) finite-difference scheme rather than the product of 
the average (P.O.A.), proved fruitful, eventhough it is generally considered 
that a non-linea'r instability can result in regions of flow reversal when 
the average of the product scheme is used. 

The problem of the treatment of boundary conditions ,at an impermeable 
wall in viscous compressible flows reduces to that of the calculation of 
the pressure or of the density. In this research the wall density was 
calculated from the continuity equation written at the wall. Peyret and 
Viviand [1 ] report that such, a technique is. of delicate use and may 'lead 
to strong oscillations or even to divergence if no artificial viscosity 
term is added to the continuity .equation; and that-, in particular, in the 
case of separated flows' negative values of the density may be obtained.. 

This, in fact, was what happened at the trailing-edge point using the 
continuity equation. This problem was overcome by using an explicit dis- 
cretization based on the leap-frog scheme, only at the trailing edge (cf. 

Eq. -4.23). 

In general it was found that at a given time-step iterative convergence 
to the tolerance of 10~® occurred in about 7 iterations. While carrying 
out these iterations the downstream boundary conditions varied with every 
iteration-. Progress in time was made by increasing the time interval At 
by .01 at the end of each time step. The problem of wiggles was overcome 
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by applying the Shuman filter at the end of every 5 time steps. Provision 
was made in the computer program to store the solution obtained at the end 
of any desired time step in a file; and the ability to read back from this 
storage file and to restart the program where it last left off was also 
incorporated. The computer program also locates and calculates the maximum 
change that occurs in a typical flow variable, such as density, along every 
i = constant line that passes through the region of the shock. Thus at the 
end of any desired time step the location of the shock and the change in 
the density across it is automatically recorded. This information is used 
again if necessary in the generation of a new coordinate system wherein the 
coordinate attraction technique is used to refine the mesh in the immediate 
vicinity of the shock. 

On an average it took 0.525 minutes of computer time (on a UNIVAC 1106) 
to achieve iterative convergence at each time step. The stand-off distance 

' - \ ' 'i' * 

of the bow-shock became quite constant after about 320 time steps each 
increment in time-step being equal to ,,01* All the programs were, however, 
run up to 400 time steps and the total computer time requirement to achieve 
this "steady-state” solution was about 2 hours and 30 minutes*, 

5*1. Discussion of Results : 

The numerical solution of the complete Navier-Stokes Equations for a 
supersonic flow: was obtained for the flow about a two-dimensional circular 
cylinder* The 'uniform flow conditions 'used in the computations -were = 4*6, 
Ke^ = 10,000, T* = 167‘'K, p* = 14,93 N/M^ and = 0.72. .For. these free- 
stream conditions the coefficient of -viscosity works out -to be - 
y* = 1.13154 X 10“^ kg/m-sec and the density p* = 3.11593 x 10“^* kg/m^..- 
The ratio of the specific heats was assumed to be y = 1.40. All calculations 
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presented in this report have been performed for the isothermal wall 
temperature T” of 556 K. The diameter of the circular cylinder was 
2R^ - 0.3048m. The Knudsen number for a perfect gas' is defined by the 
expression /Tir/2 (Il^/Re^) , which for the above free-stream conditions is 
6.821529 X 10-*^. 

The graphical results presented correspond to the steady state solutions 
at a characteristic time of 3.20, with the exception of the Mach contour 
plots which are presented at two periods of time so as to give an insight 
into the formation and progress of the bow shock wave as the solution of 
the N-S Equations advances in time towards a steady state. 

Figure 3 shows the physical field, which constituted the ^computational 
domain and Figure 4 represents the transformed field used in the numerical 
computations. A fairly compact field with 39 lines in the C^direction and 
35 lines in the n-direction was used. Even so the computer program required 
62K of core capacity on a UNIVAC 1100 series computer. 

Figures 7a-b are the Mach contour plots at the characteristic times of 
0.8, and 3.2, and depict the progression of the bow shock wave from the body 
to its steady state stand-off distance. The Mach contour interval is 0.1. 

.. i ^ M 

The sonic lines between the bow shock and the body, in- which region the flow 
is wholly subsonic, are indicated in Figure 7b. It can be seen from this 
figure that aft of the body too, a subsonic region exists which extends up 
to a distance of about 2 times the, diameter' of the cylinder from its center. 
Behind this subsonic region, the flow again is supersonic. In the field 
shown in Figure 3, the computational domain downstream of the body was limited 
by a semi-circle of radius 2,5 and the boundary conditions on this exit plane, 
as mentioned in Sect. 4.2.4, were established by solving the Euler's- Equations 
on it. Since the doxmstream boundary is located beyond the subsonic region 
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and wholly in a supersonic field of flow the use of the Euler's Equations 
is thus seen to be perfectly valid and accurate. In Figure 3 it can be 
noticed that the upstream boundary of the field starts and ends vertically 
above and below the cylinder respectively. This was dictated by the need 
of having to prescribe the free-stream conditions at least up to those 
points. 

As the n-line spacing was already sparse to begin with, the scope for 
mesh refinement in the region of the shock was. very much restricted. In 
Figure 8 the concentration of the q-lines in the region of the shock is 
exhibited. 

Figures 9 through 12 depict the variation of density, pressure, temper- 
ature and velocity from the front stagnation point to the upstream boundary 
along the symmetry line (I = 20) in the steady state. Figure 9a is the 
density distribution without mesh refinement, while 9b is the density 
variation with mesh refinement in the shock region. While the trend of 
density, pressure, temperature and velocity distributions seems satisfactory, 
the shock stand-off distance is more than what the ideal theories of Refs. 

[6] and [7] predict. This effect is first due to the solution of the Navier- 
Stokes equations and second due to the introduction of numerical dissipative 
terms. Moreover, since the present method is designated for capturing the 
shock, the shock stand-off. distance is not forced in advance as is the case 
with the shock fitting method. 

Figure 13 shows the variation of the coefficient of pressure (C ) along 

P 

* I 

the upper half of the cylinder from the front stagnation point to the 
trailing edge. 

Figures I4a-b show the distribution of wall densities from O'* to 90° 
normalized with the stagnation value without and with mesh refinement 
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respectively. Finally in Figure 15 the results plotted in Figure 14a are 
compared with the experimental results quoted in Ref. 11. It is seen that 
the numerical results of the Navier-Stokes equations , for the f ree-stream 
conditions considered, yields results which are quite in agreement with 
those obtained by experiments. 

Figures 16 and 17 show the computational domain and the sensity ratios 
for a cylinder with aft body for M =4.6 and R = 10*^. 

00 Q 
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6.0, CONCLUSIONS 


A computer code for the numerical solution of the Navier-Stokes 
equations for viscous compressible flows' with detached bow shocks ahead of 
,two- dimensional arbitrary shaped bodies has been developed. The outer up- 
stream and downstream boundaries can also be chosen arbitrarily subject 
to some obvious physical requirements. The numerical generation of 
coordinates is performed separately and then the transformed Navier-Stokes 
equations are solved on these coordinates ♦ 

The results reported in this paper pertain to circular cylinders with 
atid without an aft body as body shapes. The outer boundary is formed of 
a hyperbolic arc on the upstream side and of a circular arc on the down- 
stream side.. An implicit finite difference scheme with point-SOR is used 
to solve the Navier-Stokes equations. The computational domain is comprised 
of 39 X 35 mesh points. Iterative convergence with a tolerance of 10"^ at 
each time step of At = .01 is obtained in about 6 iterations. The total 
computer time to achieve a steady state solution on the UNIVAC 1100 series 
computer, including the time required in the generation of curvilinear 
coordinates, is about 3 hours. 

The method developed for the concentration of coordinates in the shock 
region seems to work in the sense that one or two ri“ lines are added to 
this region. However, because of the coordinate shifts, a cruder mesh 
spacing results in the non- shock region, which occasionally gives rise to 
oscillations. A better and finer coordinate spacing before the application 
of the mesh refinement is expected to yield ,better results. 
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Appendix A ' 


This Appendix summarizes the basic rules of tensor calculus [23], 

[24] used in transforming Eqs. (3.1) - (3.3), and the expanded form of the 
pertinent terms which appear in Eqs. (3.17) - (3.20). In all the formulae 
given below, repeated indices imply summation. 

Let x^ be the Cartesian coordinates and 5^ the general curvilinear 
coordinates. Then the metric coefficients are 


_ 

(A-l) 

ij _ 

(A-2) 

ik .i 

(A-3) 

g = det (g^j) 

(A-4) 

9(x^,X2,X3) 

J - /g - 

3(e^C^,C^) 

(A-5) 

The element of length ds is given by 


(ds)^ = 6 . . dx, dx, 
ij 1 J 


= g.^ dK^ de^ 

(A-6) 


Based on the above definitions, we collect other formulae from tensor 
calculus . 

Christoff el s 3 nabols: 



1 m 

2 ® 




8x 




9g 


2 ^^ 




(A-7) 
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(A-8) 


which is symmetric in j ,k. Contracting the i and j indices, we have 

k 

■■kr 28 35' 


Covariant Derivative: 


Divergence of A Vector: 


^,k ^ ggk kr ^ 


div v = ^ (v^ v^) 


Divergence of A Tensor Yielding Contravariant Components: 


/j, 1 9 , 7 - Ikv , „i rs 

(div x) = “ — (*'S "T ) + I^rs 


ae* 


Laplacian of A Scalar (j>! 


1 . a / 1 ~ ik aij> s 

= _ (/g g 


/g sr ar 

The transformation of Eqs. (3.1) - (3.3) is now direct. Using 
contravariant components, Eq. (3.1) becomes 


35' ■ 


Equation (3.2) becomes 


a . i. . 1 9 . 7- ik. . „i rs ^ 

(pv ) + — — ^ (/g T ) + X =0 


/g a? 

where i = 1,2 for two dimensions. 

Equation (3.3) becomes 

9'5' , 1 a . , i. ^ 

— — 7 (P^g b ) = 0 

' «^.35' 

The expanded form of these equations are, Eqs. (3.11) - (3.13). 


(A-9) 


(A"10) 


(A-il) 


(A-12) 


(A-13) 


(A-14) 


(A-15) 
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The terms appearing in Eqs. (3.17) - (3.20) are 


Defining 


A - i ^ 4. Z. rvl ^ \ 

® g 3n ®22 95^ g ^®22 ^^11 ^ ^12^ 


“610 (rj, + (g 


9K 


9 K v 


12 '"12 ■ * 22 '" '/- '^22 95 ®12 9n 


^ + '’ll 2 r ‘2 d *2 + r* 2 D“) 


9P 


^ g ^®12 95 ®11 9n^ *** g ^^12 ^ 22 ^ 


- g (ri + r2 )} +- (g -^ - g -^) 
^12 ''"11 42 ^^ ^ 9 ti ®12 95^ 


^ ^ 2 r |2 Dl^ + r |2 


'I' = e (*^ v^) + “ (/g v 2 ) 


. 9 , p , 9T 9T. •> 

■*■ “o 95 ^,~ ^®22 95 hi 9n^’^ 

- ' /g 


+ a 


, V /■ 9T 
^ — ^ 8 n 8 


9T 


o 9 ti ‘=’11 an ^12 95 


■)} 


the terms D'^f 


Gy - Sy/-'g 


, v^, v2, etc. , are 


2p[G22 - G^2 % ‘'' ^^22 ’’ll “ ^12 ^12^'^ 


^^2 ^12 “ ^12 ’’ 22 ^'^^ 


(A-16) 


(A-17) 


(A-18) 


(A-19) 


• (A-20) 
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(A-21) 


d12 = + G 22 Vj) - Gjj (Uj + v^) 


'*' *®22 ’’u '’' '^12 ■ '^12 ^''11 


‘S2 I2 °11 ’■22 - <=12 fil ■<■ 


““ - 2'^t'=U ’'n - <=12 '’? «=11 ■'U - °12 '■il>“ 


+ «=11 I2 - <=12 ‘■l2>’'l 


K = Muj + + (ri^ + ty-u + (rij + ry.] 


V . Ku + p(G^ G22 u + G^2 Gj2 V + u)vjj 


+ uIGj2 Gj2 u + (G,,)2v}v 2^ - y{G^, u + (G,^) 


22 ' '='12 


12^ 


- m((G22)2u + G 22 G22 v}v^2 


= Kv - 


'*'<=11 <=12 “ «12="’'>''!i 


- h{(G22)^u + Gj2 G 22 V - 


+ ^{(G^p^u + G21 G22v}v*2 


+ u {622 0^2 u + Gj2 G 22 V + v)v22 


(A-22) 

(A-23) 

2v-v>vl 

9 ^ 

(A-24) 


<A-25) 
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Appendix’ B 

Equation for Numerical Coordinate Generation 


As in Appendix A, we shall denote the Cartesian coordinates as x^ (the 
index i serving as a label only, having no tensorial significance) and the 
curvilinear coordinates either as x^, or x^ or 

It was mentioned in Section 2 that we need formulae in which' the Cartesian 
coordinates x^ are treated as dependent variables, while the curvilinear 
coordinates are treated as independent variables. To achieve this goal, 

we shall use some formulae from general tensor calculus. In all the expressions 
given below, the repeated indices imply summation. 

Formula for the Second Derivative: 

Let and be two general coordinate systems. The formula for the 
second derivative [23] is 

3^x^ _ pP 3x^ _ pt 3x^ 3x^ 

8? 8?" 3? 3i^ 87“ <B-l) 


X — X X — X 

where and" F^j^ are the Christoff el symbols in the x and x coordinate 
systems respectively. (Refer to Eq. (A-7) for definition). Since we are 
considering the transformation between a Cartesian and a curvilinear system, 

X — X 

hence either x or x is a Cartesian system; 


If x’^ is a Cartesian system, then F^^ = 0. Writing x’^ = x^ and x^ = 5^, 


we have 


3^x 3x 

- _ V-. ^ pP r 

8{“ 8CP 


(B-2) 


which is the formula for the second derivative of any Cartesian coordinate 
with respect to the curvilinear coordinates. 
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Next, if X is a 
and x^ = 5^, we have 


^ ' Cartesian system, then = 0. Writing ?■ = x^ 


111 

9x. 3x ii 3x. 3x 

Jl m £ m 


(B-3) 


which is the formula -for the second derivative of any curvilinear coordinate 
with respect to the Cartesian coordinates. 

The use of Eqs. (B-2) and (B-3) along with the equations 


9x. 

3t xp £ 


(B-4) 


„,n 

r 3g _ ^n 




(B-5) 


yields a series of useful equations. 

Inner multiplication of Eq. (E-~2) with and use of Eq* (B-5) yields 

oX< 

T,r ^ H s_ 

» f\ ' ♦ • 

_x .^3 


s 35 SC-" 


Using Eq. (B-4) in (B-7) , we have 


3x 3“^x 


s 


r5, =grt__s 

- 35*^ 35^ 35^ 

Introducing Eq* (B-4) in Eq* (B-3), we have 

= pr ip Jq j 
3x^ 3x^ ij g^p g^q 

Another form of (B-8) can be obtained by using (B-7), which is 


(B-6) 


(B-7) 


(B-8) 


3^5^ 


^ . . 3^x 3x 3x. 3x 

rt xp iq s s £ m 

= “ g g g 


35 ^ 35 ^ 35 ^ 35 ^ 35 ^ 

It must be noted that the right hand sides of Eqs. (B-7) - (B-9) have 
differentiations with respect to 5^, as desired. 


(B-9) 
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Laplacian: 


Setting A = min Eq. (B-9) , using Eqs. (A-1) and (A-3) , and 
over the index m, we get 

- - e" -I 

8| 8{ 85^ 

In two dimensions, writing 
1 

C^.= = .n, x^.= X, X2 = y 

and using 


622^^’ 


8i2/s> 


(which are a consequence of Eq. (A-3)), we have 

= [(822 y„ - 2 Si2 y^r, * ®il 


- <h2 ='« - ='cn ®11 


3/2 


where 


V n [(g22 - 2g^2 + 811 


(822 2g^2 ■'■ ®ll 


1 So<5 (Si 


3/2 


’ll ®22 '^^12 

Similarly, using (B-3) and (B~4), we easily obtain 

C ' ‘''11 *n ^ ■ ''12 % y?' + ^22 *{ 

where = C and 

The Laplacian. of a scalar function f(?^) is obtained as dlv 


72f = A. 9 


ik 3f 


V"f = (/i g--^) 




sr 


summing 

(B-10) 

(B-11) 

(B-12) 


(B-13) 

(B-14) 

(B-15) 

(B-16) 

(grad f), 
(B-17) 
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which in two dimensions has the form 


+ (2gi2 
+ (2gj_2 




a 


12 S22 ^11 


1 _ 


^ 22 ^ 


r2 

12 


^ So o ^ 


2 _ 


22 ^11 


®11 


(B-18) 


Unit Tangent and Normal Vectors: 

In the generation of coordinates we have taken clockwise traverse 
along the body contour as positive. Denoting the unit tangent and normal 
vectors as t and n respectively, and k the constant unit vector normal to the 
plane of the curve, the vectors (t, n, k) are assumed to form a right-handed 
system. The unit tangent and normal vectors for -the C = const. auH 
ri = const, curves are 



const . 






(B-19) 


(t) 


~ n = const. 




(ixg + jy^) 


11 


(n) 


1' 


5 = const. /- 
•'S- 




22 

1 


n - const. r~ 


(iy - ix ) 




11 


(B-20) 


(B-21) 

(B-22) 


where i and j are unit vectors along x and y respectively. The resolved 
parts of the velocity vector along the curves 5 = const, and n =» const, are 


(v*t)^ = 


const. 


7^ <“Sl2 '’Hz> 

’^22 


(B-23) 
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1 


(B-24) 


(v*t) 


const . 




<»Sii + 


vg^2> 


11 


(v-n)j. _ 

~ ~ 5 - 


const . 



(B-25) 


(yn) 


n 


const. 



(B-26) 


where u and v are the co'ntravariant components of the velocity vector v, 
which are related to the Cartesian components through Eq. (3.26). 
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Transformed Plane 
(Natural Coordinates) 

Figure 2. Field Transformation 










Figure 4. The Trans foraed Computational Domain. 
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Figure 7a. Mach Contours at Time 0,8. 
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Figure 8 



The Computational Domain with Mesh Refinement. 
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Figure 9a. Density Variation Across the Shock on Stagnation 
Line. 
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Figure 9b. Density Variation Across the Shock for Field with 
Mesh Refinement. 
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Figure 10. Pressure Variation Across the Shock on Stagnation 
Line. 
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Figure 11. Temperature Variation Across the Shock on Stagnation 
Line. 


52 



VELdCITT 

-UO 0.60 



Figure 12, Velocity Variation Across the Shock on Stagnation 
Line. 
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Figure 13. Coefficient of Pressure Distribution on Top of 
Cylinder . 
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Figure 14a, Density Distribution on Cylinder from 0® to 90°. 
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Figure 14b. Density Distribution on Cylinder from 0“ to 90® 
for Field with Mesh Refinement. 
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Figure 15. Comparison of Density Distribution with Experimental 
Results . 
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Figure 16. 


Computational Domain for a Cylinder with an Aft Body. 
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I'igure 17. Density Distribution on Cylinder from 0 to 90® for 
> Cylinder with an Aft Body. 
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